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Abstract 

We have studied collective modes of quasi-2D Bose-Einstein condensates with multiply-charged vor¬ 
tices using a variational approach. Two of the four collective modes considered exhibit coupling between 
the vortex dynamics and the large-scale motion of the cloud. The vortex presence causes a shift in all 
frequencies of collective modes even for the ones that do not couple dynamically with the vortex-core. 
The coupling between vortex and large-scale collective excitations can induce the multi-charged vortex 
to decay into singly-charged vortices with the quadrupole mode being one possible channel for such 
a decay. Therefore a thorough study was done about the possibility to prevent the vortex decay by 
applying a Gaussian potential with its width proportional to the vortex-core radius and varying its 
height. In such way, we created a stability diagram of height versus interaction strength which has 
stable regions due the static Gaussian potential. Furthermore, by using a sinusoidal time-modulation 
around the average height of the Gaussian potential, we have obtained a diagram for the parametric 
resonance which can prevent the vortex decay in regions where static potential can not. 
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I. INTRODUCTION 


The dynamics of a trapped Bose-Einstein condensate (BEC) containing a vortex line at 
its center has been the object of our studies. We have studied the effects of a multi-charged 
vortex in free expansion dynamics. These central vortices contribute with the quantum pressure 
(kinetic energy) which increases the expansion velocity of the condensate |23| . Consequently, 
our work culminates in describing the collective excitations of a vortex state as well. Here the 
vortex-core dynamics couples with the well known collective modes |22| . Furthermore, we shows 
that it is possible to excite these modes using modulation of the s-wave scattering length. Such 
a technique has been already applied to excite the lowest-lying quadrupole mode in a lithium 
experiment [19]. The motivation for these works is the possibility of experimental realization. 
Now our focus is in the anisotropic oscillations of the vortex-core. In other words, oscillations 
that lead the vortex shape from circular to elliptical. Such deformation is a symmetry breaking 
of vortex state, and can result in changes of dynamical stability. 

The presence of vortices in condensates can also shift the frequency of collective excitations. 
The frequency shift of quadrupole oscillations have been analytically explored for positive scat¬ 
tering lengths by using the sum-role approach BH. as well as the effects of lower-dimensional 
geometry in the frequency splitting of quadrupole oscillations [Tj. 

First of all, multi-charged vortices in trapped ultracold Bose gases are thermodynamically 
unstable, which means that a single Acharged vortex tends to decay into £ singly quantized 
vortices. Thus the configuration of separated singly-charged vortices has lower energy instead 
a single vortex with the same angular momentum. Although such a state with multiple singly- 
charged vortices is also thermodynamical unstable when compared with a vortex-free conden¬ 
sate. These multiple vortices spiral outward from the condensate until remain only the ground 
state. 

The vortex dynamic instability has so far been studied in the context of Bogoliubov ex¬ 
citations mmm- Indeed, the vortex state possess certain Bogoliubov eigenmodes which 
grow exponentially and become unstable against infinitesimal perturbations |8]. These vortices 
present several unstable modes being a quadrupole mode the most unstable. For instance, let 
us consider the work in Ref.jHJ. There the authors studied the modes of quadruple-charged 
vortex. Among of them, only three modes are unstable. These unstable modes have complex 
eigenfrequencies (CE) and are associated with /-fold symmetries. These symmetries are: 

• Two-fold symmetry; the quadruply-charged vortex splits into four single vortices arranged 
in a straight line configuration. 
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• Three-fold symmetry; the quadruply-charged vortex splits into four single vortices ar¬ 
ranged in a triangular configuration, i.e., there are three vortices forming a triangle with 
each vortex representing a vertex. The fourth one is at center. 


• Four-fold symmetry; the quadruply-charged vortex splits into four single vortices arranged 
in a square configuration with each vortex placed in a vertex. 

Our target is to describe them as a result of the coupling between the vortex-core dynam¬ 
ics with the collective modes of the condensate. In order to achieve this goal we have used 
variational calculations focusing on the description of only one of the unstable mode (specially 
two-fold symmetry). The variational description becomes very complicated as we increase the 
number of parameters. Fortunately the most relevant unstable mode is also the easiest one to 
calculate within the variational approximation. 

Furthermore, there are some works which add a static Gaussian potential centered in the 
core of a vortex with a large circulation which results into a stable configuration for the multi- 
charged vortex H2H3]. Based on these works we checked the dynamical stability for a static 
as well as dynamic potential due to a Gaussian laser beam placed in the vortex-core, when 
compared with the multiple vortices state. 

This paper is organized as follows: In section [TT| the quasi-2D approach is introduced. We 
discussed the wave-function used with the variational method in section IIIII and detailed the 
calculation of the Lagrangian in section IV Section [V] contains equations of the motion and 
their solutions, i.e. the stationary solution, collective modes, and the fully numerical calculation 


of Gross-Pitaevskii equation (GPE). In section VI we made a dynamical stability diagram 


considering a static Gaussian potential while in section VII we made a parametric resonance 
diagram due to a dynamical Gaussian potential where its height is sinusoidally time dependent. 


II. QUASI-2D CONDENSATE 


The presence of a large number of atoms in the ground state allows us to use a classical 
field description [,18j. Where the non-uniform Bose gas of atomic mass m and s-wave scattering 
length a s . The scattering length is smaller than the average inter-particle distance at absolute 
zero temperature. Its dynamics is given by the Gross-Pitaevskii equation m 


dt 


~V 2 + V(r) + I/ 0 |T(r,f)r 
2m 


T (r,t), 


where the interaction strength between two atoms is 

47 rh 2 a* 


Vn = 


m 


( 1 ) 


( 2 ) 
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In order to suppress possible effects due to motions along the axial direction, we consider a 
highly anisotropic harmonic confinement of the form 


= \rnu>y +^muj 2 z z 2 , 


V (r) = V ± (r ± ) + V z (z) 

( 3 ) 

with l o z c o p . With this condition the condensate wave-function can be separated as a product 
of radial and axial functions, which are entirely independent. This yields a quasi-2D Bose- 
Einstein condensate |H EDI ESI, and leads to 

V(r,t) = N$(r ±1 t)W(z,t), (4) 

where 


W(z,t) = 


d z \J 7T 


exp 


iu z t 


2 dl 


By replacing (J3]) and ([f]) with (J5]) into the Gross-Pitaevskii equation ([Tj) , we obtain 




mw m + 2 


= 


h 2 

■—Vi + 


h 2 


h 2 z 2 


+ V (r) + NU 0 |$W| 5 


<f>W. 


(5) 


( 6 ) 


2m ' " L ' 2 mdl 2 mdj 

The product $ (r_|_,£) W ( z,t ) is normalized to unity, thus the number of atoms appears multi¬ 
plying the coupling constant Uq. Now we can multiply Eq. ([6]) by W* ( z,t ) and integrate this 
equation over the entire z domain. Since h 2 /2md 2 = Hou z /2, and h 2 /2md^ = muj 2 / 2 we obtain 
the following simplified equation 


<9$ (r L,t) 


+ i4 ( r± ) + m j 2D | 4 (r 

2m 




where 


Uon — 


Un 


= 2y/2n 


h 2 a , 


(7) 


( 8 ) 


d z y/2n ” " md._ 

Let us then write the Lagrangian density which leads to quasi-2D Gross-Pitaevskii equation 
([r]) for a complex held $ (rj _,t) normalized to unity. So the Lagrangian is given by 


r _ id 

L-in — —— 


2 

+ ^ | Vi<h (r_L, t) | 2 + V ± (r_i_) |$ (rj_, t)| 2 + |<f> (r ± , t)\ 4 . 


T i , .d$(r±,t) . . <9$* (r_L, t) 


dt 


dt 


(9) 


III. BREAKING WAVE-FUNCTION SYMMETRY 


In order to examine the coupling between the vortex-core dynamics and the collective modes 
as well as their stability, we choose the situation where a multi-charged vortex is created at the 
center of a condensate. Its wave-function can be written in cartesian coordinates as 


(rj_, t) oc 


/i i r/r 

X 

2 

y 

l [x/ix(t)f + 2xy/£ xy (t) + [y/iy(t)} 2 + lj V 

R x (t)_ 


Ry (t) _ 


D iS(r ± ,t) 


( 10 ) 


4 




























The sizes in each direction are given by R { (t). They are known as Thomas-Fermi radii, since 
the wave-function vanishes for x > R x and y > R y . The vortex-core sizes are given by (t). 
They are of the order of the healing length for a singly charged vortex. The parameter £ xy (t) 
is responsible for a complete description of the quadrupole symmetries between vortex-core 
and condensate. The wave-function phase S (r, t) must be carefully chosen within the context 
of the variational method. Because the phase must contain the same number of degrees of 
freedom as the wave-function amplitude. Since we have one pair of variational parameters for 
each direction in the wave-function amplitude (£ ?; and Ri ), we also need a pair of variational 
parameter in the wave-function phase (Bi and (7*): 

2 2 4 4 

S (r ± , t) = £ arctan 0) + B x (t) y + B xy (t) xy + B y (t) y + C x (t) y + C y (t) y. (11) 


Thus Bi ( t ) and Ci (t) compose the variations of the condensate velocity field allowing the 
components (t) and Ri (t) to oscillate with opposite directions. While B xy (t) gives us the 
contribution of the distortion £ xy ( t ) for velocity field which changes the axis of the quadrupole 
oscillation. Note that we are not using a parameter which yields a scissor motion to the external 
components of the condensate, since it has already been shown that such a motion is not coupled 
with neither breathing nor quadrupole modes mm. 

This choice for our wave-function implies that our vortex-core might have an elliptical shape. 
It is enough to destabilize a multi-charged vortex and allow it to decay splitting itself into several 
vortices, each one with unitary charge. 

Following the variational method used in Ref. [I5l [T6l 22J [23], we substituted (10) into 
§. and performed the integration over the spacial coordinates, L 2 d = f £- 2 Although 
the Lagrangian density (J9]) cannot be analytically integrated since it does not keep the polar 
symmetry. One way to proceed is to introduce small fluctuations around the polar-symmetry 
solutions into the wave-function, and to then to make a Taylor expansion. Thus we can take 
advantage of the approximate polar symmetry of the vortex-core while the fluctuations act 
breaking the vortex-core symmetry. These calculations are discussed in detail in the next 
section. 


IV. EXPANDING THE LAGRANGIAN AROUND THE POLAR-SYMMETRY SO¬ 
LUTION 


Within the Thomas-Fermi approximation the trapping potential shape determines the con¬ 


densate dimensions. The wavefunction (10) is approximated by an inverted parabola except for 
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the central vortex. So that its integration domain is defined by 1 — x 2 / R 2 — y 2 / R 2 > 0. Some 
care should be taken when calculating the kinetic energy |Vj_<f> (r_|_,t)| 2 before integrating. The 
vortex presence inserts an important term in the gradient, while the rest of the gradient is 
neglected in the Thomas-Fermi approximation. That means that the density varies smoothly 
along the condensate except in the vortex. 

By introducing deviations from the equilibrium position in our parameters 

(t) ~ £o + Stj ( t ), (12) 

Rj (t) Rj R 0 + 5Rj ( t ), (13) 

we can expand in Taylor series the deviations of the Lagrangian. In this way we have 

L = L (0) +L (1) + L (2) + ... (14) 


The linear approximation is obtained by truncating the series in second order terms, this leads 
to: 

• Terms of zeroth order in L ^ being responsible for the equilibrium energy per number of 
atoms. 


• Terms of first order L'- l > that vanish due to the stationary solution of Euler-Lagrange 
equations. The equilibrium configuration has polar symmetry. 

• Terms of second order L^ carries the collective excitations. Their Euler-Lagrange equa¬ 
tions result in a eigensystem whose eigenvectors are composed of deviations (SRj and 

Notice that B j and C, from phase ([IT]) also must be considered as first order terms since they 


lead to deviations in the velocity held. In order to evaluate all the necessary integrals in Eq.(14), 
it is convenient to use R (£) /R t (t) = cc* (t) instead of (t). This change is explained due to all 
these integrals result in functions of cto = £o/Ro- Thus, a we use a* (t) ~ «o + Sai (t) instead of 
6 ('t ) « £o+6£i ('t ). The same happens for a xy (t) = R x (t) R y (t) /£ xy (t) where a xy (t) Rj 5a xy (t). 
Hereafter we omit the time dependences for simplicity, and we named zeroth order functions 
as Ai = Ai (£,a 0 ) as well as the other integrated results as R = R (£,a 0 ). Such functions are 
described in Appendix [Aj 

The proportionality constant in wave-function ( |Toj ) is found through normalization, being 
N 0 = R x 1 R y 1 [Aq + R {5ot x + Scy y ) + R (fioi 2 + Sa 2 ^j + R5a x Soiy + J^cr 2 ^] 


« o 1 




(Sa x + 8a y ) + ^ ~2 ~ ($ a l + d a l) 


_ h\ A c _ Aa 2 
A 2 A 0 ) SaxSay A 0 Soixy 


(15) 
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By calculating the Lagrangian integrals we obtain 


p 2 |$| 2 dr ± = N 0 R x R y < R x 


R: 


—i 


A\ + I§ ( 8 a x + — 8 a y ) + + I 7 8 a 2 + I% 8 a x 8 ay + Ig8a 


A\ + R ( — 8 a x + 8 a y ) + I 7 8 a 2 x + R 8 a 2 + I% 8 a x 8 a v + Ig8a 


d>*-<|>- 

dt dt 


dr 1 — N 0 R x R y R X B X Ai + J 5 ( 8 a x + — 8 a y 


Y xy 


A2B xy I w 8a xy + R^yBy 


A\ + I 5 ( - 8 ai x + 8 a y 


- 2 KC x 

-\KCv 


A 2 ~b in ( 8 oi x -b —801 


5 ~ v 


Ag + -in I ~^8a x + 8a y 


~xy 


(16) 


(17) 


I'V±®| 2 dr ± = N 0 R x Ry {A,R^ {11= + 2B% + Ip) + 27l 2 fl* I BJX + B„C y ) + -4 3 B® 


and 


i 2 

A —j [A4 + Ii 2 8a x + IizSoty + 1148a 2 + Ii^8a 2 + I\g8a x 8a y + In8a 2 y ^ 

i 2 

+ -— [A4. + Iiz8oi x + I\28a y + Ii^Sa 2 + Jii^o;^ + I\g8a x 8a y + In8a xy j 

Ky 


+ 


£ 2 

Rl L 


A 5 —— ( 8 R X + 5i2y) + ■—2 (<5ii^ + £ifc + 8 R x 8 R y ) 

-0-0 -fin 


L o 


| 8 R x 8 a x ~b 8 R x 8 a y 8 R y 8 a x 8 R y 8 a y 

rig \ o o 


— 77-^18 (^®a: + 5o(y) + /19 + <5o^) + I2o8oi x 8a y + /2i<5<a 


(18) 


I |$| drj_ — NqR x R,y [Ag + I 22 ( 8a x + + I 23 ( 8a 2 + 8a + I2i8a x 8a y + . (19) 

By scaling according to Table [TJ each of the three first terms from (14) are given by 

V2npA 6 


L<°> = An 1 


£ 2 ( 1 

Al^ n n H- TT ( A 4 + -Ac 


' po 


r pO A 0 


( 20 ) 


= -A 


-1 


A i r P o yBx + ByJ + ~A 2 r p0 [C x + C y 


~\~Sp ( SR X + 5Ry ) “f" aSq, ? 


L( 2 ) = -An 1 


• ( St \ • /" St 

AiTp 0 f3 x I 2-b Fi8a x + F 2 8a y J + Air 2 0 j3 y I 2—- + F 2 8a x + F\8a l 

\ r p 0 J \ r pO 


1 • ( St \ 1 • ( St 

A~A 2 r p0 ^ x I 4-b F :i 8 a x + F A 8 a y I + -A 2 r p0 ( y I 4—- + F 48 a x + F 3 8 a 1 

z \ r pO J 1 \ r P o 


+Air 2 p0 (Pi + Py) + 2 ^ 2 ^p 0 (/^xCx + ^j/Cj/) + A 3 r 6 p0 (( 2 + ( 2 ) + f/ p + <5r^) 
T Vpp8r x 8r y ~b Vq, ( 8a x T 8a y ^ T ~b Vp a ( 8r x 8a x ~b 8r y 8a y ) 


+hap ( \8v x 8oi y ~b 8v y 8a x ) ~b 2 ^po ( Ai/3 ~b -iio (3xy8a xy j T V xy 8a 


( 21 ) 


( 22 ) 
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where 


Q-9/1 ^ (0 4 A ) 2 V / 2tT7^6 

S p — 2 A 4 r p0 — - 3 - ( 2 A + A) — 


po 


r 3 p 0 A 0 


S a — A' r po (A + A) + ~ 2 ~ [A (F 5 + Fq) — A (A + F$)] + 

r pO 


2v / 27T7^6 rp 

r 2 4 n A, 

r p0 /i 0 


v = 4 1 + il(3A 1 + 4 5 )-' 


pO 


Ap - IT A 
1 pO 


AA 


2\/27r7A 

r pcA 


( 23 ) 

( 24 ) 

( 25 ) 

( 26 ) 


K = 


/fi +17 


V n = 


^I'pO ^ 

A 1 

— z 

h 

to 

1 14 + 7 i 5 

"s> to 
0 

L a 

. A^ 

-A 

A 

1 r 2 
' pO 

.A 

A 

2 v / 27 T 7 A 

I 23 

r 2 P o A o 

.A 

A i r % ~ 

^6 + h r 

A 1 

z 

A 4 £ 2 

_ 1 _ 

1 14 + -A 

"i 

"£> to 
0 

L A 

, AA 

I 19 

_ A 

-\ 

■i 

-£> to 
0 

_ A 5 

A 

2 v / 27 T 7 A 6 

A 

1 r pO A 0 

.A 

2 Ar-p 0 

A 
.A " 

A 

A 

A 4 e 

+ 2 r 2 
' p 0 

r a 

.A 

J s 

, AA 

120 

A 

"£> to 
0 

.A 

A 

2 v / 27 T 7 A 

A 

r p 0 A 

.A 


- 2 I - A 5 + 

-7 ( f 7 + p 8 

-2^--^- (2F 9 + ^~ 

A A V A 


A _ii 
A A 


- 2 £ - 7 (fi + *> 


r (Fj+ 

24-0 

>A_ii 
A A 


A 

-An 


h ( 

A ( 

_ A 
io A 

A ( 

A 

h 

2 — - 

A 


,A 

A 


— 2 A\r poA 3 - (2 A A — A A) 


pO 


V ap — 2 A'PpoA — 3 - (2 A A — A A) — 


po 


2F9+ t 

2\Z27T'jA 6 
r%A 0 

2\[2t\^Aq 

r p 0 A 


9, 


9, 


( 27 ) 


( 28 ) 


( 29 ) 

( 30 ) 

( 31 ) 



V xy = 2A\r 


Ad 2 


P o 


with 


J 9 I A \ ~Ad 2 


f-f 1+2 
Ai A 0 


hi 

A-5 


h 

Aq 


' pO 

2\p2d^Ao 


h 7 

A 4 


Fi = 


F, = 


F* = — - — 


F d = 


F q = 


T 

r 2 

P O' 

h 

h 

A\ 

Ad 

h 

h 

3Ai 

Ao 

hi 

h 

A 2 

V 

In 

h 

5 A 2 

Ao 

I 12 

h 

a 4 

Ad 

h 3 

h 

a 4 

Ad 

ha 

h 

A 5 

Ad 

ha 

h 

3A 5 

Ao 

1 22 

2 Jl 

Aq 


7 r p 0 

are d 


F* = — - — 


Fp = 


F 7 = — - — 


F« = 


h 

Aq 


— - 2 — 
Ar An 


(32) 

(33) 


(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 


the multi-charged vortex, and the interaction parameter in dimensionless units is given by 
7 = Na s /d z . In the next section, we discuss the Euler-Lagrange equations for the deviations 
that lead to the four collective modes. Where one of them is dynamically unstable. 


V. ENERGY PER ATOMS, COLLECTIVE MODES, AND INSTABILITY OF A 
QUADRUPOLE MODE 


First, in order to calculate both energy per atoms and collective modes, we need to know 
the equilibrium points r p o and «o. They are obtained from Euler-Lagrange equations for Sr j 
and <5cq, resulting in 


S p = 0, and S a = 0. (43) 

Thus we have different pairs of r p o and «o f° r each value of i and 7 , which are obtained by 


applying Newton’s method to solve these coupled stationary equations (43). Note that for 7 = 0 
its solution is trivial, given by 

r p0 = 2(2 /tt) 1 /8 7 1/4 . (44) 
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Dimensionless scale 

t 

uj-H 

h 

kujpjl 

n 

UpQ 

Ro 

dpv po 

SRj 

d p 5rj 

£0 

d P r^ 0 


d p 5r u 

B, 

dfPj 

Cj 

d~ P % 

zu 

UJ p CaJ 


Dimensionless parameters 
7 = Na s /d z 
a = £/R 


Table I: Scale table. 


These equations (43) do not have physically consistent solutions for low values of 7 depending 
on the value of £, as can be seen in fig|TJ We have evaluated the values of the pair r p 0 and 
«o for the vortex-states with l = 2,4, 7, where the lowest values of interaction are around 
7 = Na s /d z = 29, 76,125, respectively. 

The energy per atom Zd 0) increases proportionally to 7 1 / 2 being more evident for the vortex- 
free state {£ — 0 ), where 

L (0) = 4(2/tt) 1 / 4 7 1 / 2 /3. (45) 


We show this behavior for others values of £ in figj2j The energy gap between the vortex-free 
state and the remaining states corresponds to the amount of energy needed to create the l- 
charged vortex states. For instance, if a focused laser beam is used to stir a Bose-Einstein 
condensate in order to nucleate vortices, the stirring frequency must exceeds a critical value [? 
], which is defined by difference of energy between the vortex-free state and the singly vortex 
state. 

Calculating the Euler-Lagrange equations from we obtain ten coupled equations, being 
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(a) i = 2 



(b) l = 4 



(c) l = 7 


Figure 1: (Color online) Equilibrium point of parameters (r p o, and «o) by atomic interaction. 
Solid (black) line represents r p0 , and dashed (blue) line represents chq. Both are calculated 


from (43) where ao must be smaller than r p o and near to zero value. This approach shows 
itself valid for Na s /d z > 29 ( Na s /d z > 66, and Na s /d z > 125) when we have i = 2 {t = 4 and 

e = 7). 
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N a s / d z 

Figure 2: (Color online) Energy per atom as function of interaction parameter. 


five equations for phase 


5r x 

r p o 

K 

r P o 

8r x 

r p o 

St 
U l y 

r p o 


+ 

+ 

+ 

+ 


2 

F\ 

2 

Fs 

4 

F\ 

4 


<5<Te + 
<5«2; + 
<5<Te + 
8a x + 


2 

^Sa v 


Fa,- 

4 ° a y 

f 3 ,- 

~j 6a y 


I io 5cx X y 




Py + ~ 


A 2 


B + — 

Px+ a 2 


fiy+lT 


A 3 
^2 
2^4i (3 xy , 


’’pot, 

r />()Cy: 

rjloCz, 

r poCpi 


(46) 

(47) 

(48) 

(49) 

(50) 


and other hve equations for variational parameter in the amplitude 

AiT p0 /3 x + Azt^qC'X + 2 V p 5r x + V pp Sr y + V pa 5a x + W p <5ay=0, (51) 

A±r p0 /3 y + A 2 T 3 p0 ( y + V pp Sr x + 2V p Sr y + V ap Sa x + (52) 

741^(4^1+4^) + 2 (4-F3 + Cy4) +4(oo^ r a; + K t p5r y +2V^5Q; iT + K,,Q,5ay=0, (53) 

t4 1 r po (4 F‘i +4 -Fi) +-^r^o ( 4 F 4 +C 2 /F 3 ) +I^p5r x +V), Q 5r y +V^ a 5o( x +2V r Q (5ay==0, (54) 

r 2 p0 ho$xy + V^tfa-BjpO. (55) 


We can reduce these ten equations into 4 coupled equations plus one uncoupled equation. The 
equation for Sa xy is uncoupled from the others according to 


SOi X y T 


‘ZAlVxy 

P r 2 
J 10'pO 


SOixy 0, 


(56) 
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i.e., the motion represented by the deviation Sa xy is independent of the other collective modes. 
Those four equations lead to the linearized matrix equation 


MS + VS = 0, 


/ M p 0 M pa M ap \ 


^ Sr 


(2 Vp Vpp Vpa Vap\ 


( Sr 

0 M p M ap M pa 


Sr y 

+ 

Vpp 2 Vp Vap Vpa 


St 

Ul y 

M pa M ap M a M aa 


Sa x 


Vpa V ap 2V a V aa 


Sa x 

\M ap M pa M aa M a ) 


\SdiyJ 


\V ap Vpa Vaa 2 V a J 


\Sa y J 


where the entries in the matrix M are given by 

M p 
M a 

M aa 
M pa 
M ap 


= 2A U 

^1^2 r p0 


2 (^2 — AiA 3 ) 

^ 1 ^ 2 r p 0 

= 2 (A 2 - A^) l 
— AiFyr p0 , 


F\F 3 + F2F4 
F1F4 + F2F3 


F 2 
r 3 


F 2 

r 4 


4 

F3F4 


4 + ^ 

2AlA 3 


A 2 


-F\F 2 


— Ai F 2 r p0 . 


(57) 


(58) 

(59) 

(60) 
(61) 
(62) 


Matrix V results from the energy part of the Lagrangian, i.e. from Eqs. (16), (18), and (19). 
This determinant may be either positive or negative reflecting the system stability. In the other 
hand, the determinant of M cannot be negative or zero, since it results from our choice for 


the wave-function phase. The equation (57) seems the Newton’s equation therefore we can say 
that matrix M has an effect of mass-like, and matrix V works as a potential [9f. Solving the 
characteristic equation, 

det (. M~ l V - w 2 I ) = 0, (63) 


results in the frequencies of the collective modes of oscillation. Eq.(63) is a quartic equation of 
w 2 . This means that we have four pairs of frequencies being one pair for each oscillatory 
mode. Among these four modes, two of them have a static vortex representing the collective 
modes for cloud: they are the breathing mode B c , and the quadrupole mode Q c . In other words, 
these modes are similar to collective oscillations of the vortex-free state, where the difference is 
in a small shift in their frequencies depending on the charge of the vortex, as it is shown in Figj3j 
Therefore B c decreases the frequency value while Q c has the opposite effect shifting to higher 


frequency value. Note that for a vortex-free condensate £ = 0, Eq.(63) is a quadratic equation 
in w 2 . That means the system presents only two modes ( B c with w — 2, and Q c with w = \fF) 
in absence of vortex, whose frequencies are constant with respect to the interaction parameter 
7 . There are still other two modes which couple vortex dynamics with collective modes. They 
are another breathing mode B v and another quadrupole mode Q v . In this breathing mode B v , 
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the vortex-core sizes oscillate out of phase with cloud radii, while Sr x (8 a x ) and 5r y (8 a y ) are in 
phase. In the quadrupole mode Q v , both these sizes <5cq and 5Ri are oscillating in phase while 
5r x (Sa x ) and 5r y (8a y ) have a 7r-phase difference between their oscillations. These modes are 
sketched in figj4j The second quadrupole mode Q v has an imaginary frequency (fig.4c), i.e. Q v - 
mode is one possible channel to a multi-charged vortex decay into unitary vortices. Therefore, 
the multi-charged vortex decay can be explained by the appearance and growth of this unstable 
quadrupole mode due to quantum or thermal fluctuations. These fluctuations work inducing 
collective modes, which are coupled to the vortex dynamics through their sound waves. 

This model is completely consistent with CE Bogoliubov modes for l = 2, which are com¬ 
posed by only the CE mode associated to two-fold symmetry being our quadrupole mode Q v 
[ IT] . However when l > 2 this calculation is incomplete since we considered only breathing and 
quadrupole modes. Hence for a complete description it is necessary to add others symmetries 
for each higher order of £, which is not a trivial task. Because the Ansatz requires more degrees 
of freedom, that means we should increase the number of variational parameters. 

In order to check our results we proceed the full numerical calculation of the Gross-Pitaevskii 
equation (with the usual phenomenological dissipation e used since Ref. |24j). The reason of this 
dissipative description is the prevention of non-physical waves created by the grid edge. The 
initial state is calculated by evolving a trial function in imaginary-time with the parameters 


given by the equilibrium point from Eq. (43). We introduce the eigenvector from Eq. (57) 
corresponding to the unstable quadrupole mode ( Q v ). This trial function is given by 

[ x / (£o + <5£x)] + i [; y/ (£o + <5£y)] 


oc 


\l[ x l (£o + + [y/ (£o + <5£?/)] 2 + l 


r 

X 

2 

V 

/ » 

_(Ro + SRx). 


. {Ro + &Ry) _ 


(64) 

Furthermore we have done the evolution in real-time where we could check the multi-charged 
vortex decaying to an initial state containing only the deviations of (5,,-mode. In figure [5] is 
shown the evolution of the condensate in real-time for a doubly-charged vortex, such that it 
starts to split around u p t = 20.2. In figure [6] we notice that the life-time of quadruply-charged 
vortex is around u> p t = 22. It is necessary to observe that these life-times are different depending 
on the amplitude of deviations and imaginary-time evolution. It is also possible to induce the 
decaying by shaping an anisotropic trap, however our semi-analytic approach is valid only for 
an isotropic trap. 

It is interesting to observe the way in which multi-charged vortices decay by Q„-mode excita¬ 
tions, which makes the multi-charged vortices split into a straight line of vortices with unitary 
angular momentum. For instance, we see in figure [6] the quadruply-charged vortex splitting 
into four vortices and forming a straight line, then evolving based on its interaction with the 
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Figure 3: (Color online) Frequency as functions of interaction parameter with respect to 
cloud’s collective modes in (a) and (c). These two modes have real frequencies in domain of 
positive interaction ( Na s /d z > 0). Schematic representation of each collective mode is in (b) 

and (d). 


velocity fields until the final configuration. 


VI. STABILITY DIAGRAM DUE TO A STATIC GAUSSIAN POTENTIAL 


Some articles on numerical simulations propose to stabilize an multi-charged vortex by turn¬ 
ing on a Gaussian laser beam at the middle of the vortex-core. It means basically that we need 
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(c) Quadrupole mode ( Q v ) (d) Quadrupole 

mode (Q v ) 


Figure 4: (Color online) Frequency as function of interaction parameter for collective modes 
coupling the dynamics of the vortex-core with the oscillation of atomic cloud radii. Only the 
quadrupole mode (c) is unstable with imaginary frequency. Schematic representation of 
collective modes are shown in (b) and (d). B v mode has the vortex core oscillating out of 
phase with cloud radii. Q v mode is a quadrupole oscillation where vortex core is in phase with 

cloud radii. 


to add an external potential with Gaussian shape to the harmonic potential, i.e. 


Vl (rj Vtrap (r±) + V G (r_i_) 

= \ m u 2 p P 2 + ^Vbe -pa/ *o, (65) 

where the Gaussian width must be proportional to the vortex-core radius (w = V2£o)- An 
apparent objection to our approach could lie on the fact that optical resolution limit of a laser 
beam is around of some microns, while single-charged vortex core is usually smaller than 0.5/ini. 
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Figure 5: (Color online) Time evolution of the density (a,b,c) and phase (d, e, f) of 
condensate with a doubly-charged vortex. We have used jx = 20.198, Na s /d z = 100, 
e = 0.001, and a factor of 0.01 multiplying of the amplitude of deviations. 


However, multi-charged vortices may attain much larger sizes depending on its charge, the trap 
anisotropy, number of atoms, and atomic species. For instance, a quadruply-charged vortex 
in a 85 Rb condensate (N = 10 5 , a s = 100a 0 , c o p = 10Hz, and c o z = 100Hz) has 5.9 n r m. By 
applying a Gaussian beam with w = 10 fim inside of this vortex, its radius grows to 7.1pm. 
Thus we use this procedure in our semi-analytical method in order to draw a stability diagram, 














Figure 6: (Color online) Time evolution of the density and phase of condensate with a 
quadruply-charged vortex. We have used fi = 45.9552, Na s /d z = 520, e = 0.001, and a factor 

of 0.001 multiplying the amplitude of deviations. 


and show that it is enough to stabilize the quadrupole mode Q v . So we have to calculate now 
the Lagrangian part corresponding to the Gaussian potential, 


Lg — 


Vd(r±) |$(r_L)| J dr±- 


( 66 ) 
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By expanding in Taylor series the integral of Gaussian potential Vjxj (rj_) we have 

J e- p2 /% ($*$) dr ± = N 0 R x R y 

T ho { 8 ot x + Socy) + ho ( 8 a 2 T $ a y) + h\ 8 a x 8 a y 

T (8 R x 8 a x T SRySoty) —— (8 R x 8 a y T 8 R y 8 ah) T 1348 a 

Ro ' Ro 


A 7 + ^ ( 8 R x + hi?,) + (<$i£ + + ^SRJRy 

Kq n 0 


Rl 


( 67 ) 


where Lagrangian part becomes 

L g = 


Vh^7 f 1 + _^26_ + 5R ^ + ^27 
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Notice that we have terms of first order in deviations in Eq. (68), it means that the stationary 


solution is modified when the condensate is under the influence of a Gaussian potential. The 


first-order contribution in (68) becomes 


Sc* 
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(69) 


while the second-order terms are 
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(70) 
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Figure 7: (Color online) Diagram of magnitude of pinning potential by atomic interaction for 
vortex with £ = 2 and £ = 4. Hatched region represents stable eigensystem meaning the 

vortex-core become stable. 


The equilibrium points are changed to 


S P + s p = 0, and S a + .s a — 0. 


(71) 


Each terms in (70) adds a contribution to a different element in the matrix V of the linearized 


Euler-Lagrange equation (57) which then becomes 


M5 + (V + V G ) 5 = 0, 


(72) 


where 


^G = 


VnA 


0^7 

2 An 


(2 p p Ppp PpOL Pap ^ 
Ppp 2 Pp Pap Ppa 
Ppa Pap 2 Pa Paa 


(73) 


\Pap Ppa Paa 2p a J 

Since the stability of the eigensystem depends only on the Q v - frequency, we can build a 
stability diagram of Vo /hu p versus Na s /d z . In figjTj this diagram is shown considering two 
cases, £ = 2 and £ = 4. As the angular momentum £ gets larger the stable region decreases. 
Hence the pinning potential can prevent the vortices from splitting for some values of Vo /hco p 
depending on Na s /d p . 

In order to validate these stability diagrams, we make a numerical simulation of the Gross- 
Pitaevskii equation. When the Gaussian potential is turned on, we have seen that it provokes 
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some phonon-waves on the condensate surface and increases a little the vortex-core size be¬ 
sides preventing the vortex decay. Figure [8] shows phonon-waves rising and vanishing due to 
dissipation. The same phenomena may be seen in figure [9j 

The vortex decay happens when the sound waves couple the quadrupole mode from the edge 
of the condensate with the vortex-core, which breaks the polar symmetry of vortex. Therefore, 
the pinning potential acts as a wall reflecting these sound waves, and preventing the vortex 
symmetry break. 


VII. DIAGRAM OF STABILITY DUE TO A DYNAMIC GAUSSIAN POTENTIAL 


In section VI, we have seen that it is possible to make a multi-charged vortex stable using 
a static Gaussian potential. In addition, we calculated a diagram of height versus interaction 
strength which shows the stable region. Here we propose to stabilize a multi-charged vortex 
with a sinusoidal modulation of height of the Gaussian potential with an amplitude given by 

SV, 

Vo (t) = V 0 -SV cos (Of), (74) 


at the specific region of interaction strength where the static potential is not capable of stabi¬ 
lizing the vortex, i.e. 0 < Na s /d p < 160. The equation for this case is given by 


MS+<V + V G 


1 — — cos ( Of 
Vo 


6 = 0 , 


(75) 


where matrices M and V can be found at Eq.(57), and Vlg is given by Eq. (73). By scaling the 
time as follows 

n~ 

—t -a 2r, (76) 

UJp 

we obtain the Mathieu equation 


A 2 C ct r 'x 

—<5 + | A - 2—Qcos(2r) U = 0, (77) 

where A = M -1 (V + Vlg) and Q = (1/2) M~ 1 Vlg are constants depending on the initial 
conditions. This equation becomes solvable by using Floquet theory p El EDI EH ES]- The 
basic idea of this theory is that if a linear differential equation has periodic coefficients, the 
solutions will be a linear periodic combination of functions times exponentially increasing (or 
decreasing) functions. Thus linear independent solutions of the Mathieu equation for any pair 
of A and B can be expressed as 

^ ( r ) = e ±VT P (±r), (78) 
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Figure 8: (Color online) Time evolution of the density (a, c, e) and phase (b, d, f) of 
condensate with a doubly-charged vortex. We used id/hjjp = 20.198, Na s /d z = 100, 
Vo/fadp = 150, e = 0.001, and a factor of 0.01 multiplying the amplitude of deviations. 

where rj is called the characteristic exponent which is a constant depending on both A and Q, 
and P (r) is 7r-periodic in r that which can be written as an infinity series 

OO 

S(r) = e^ 6 2ne 2mT , (79) 
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(f) oo p t = 100 


Figure 9: (Color online) Time evolution of the density (a, c, e) and phase (b, d, f) of 
condensate with a quadruply-charged vortex. We used /i/hiUp = 45.9552, Na s /d z = 520, 
V 0 /huj p = 500, e = 0.001, and a factor of 0.001 multiplying the amplitude of deviations. 


with b 2n being a Fourier component. Doing the substitution of (79) into (77), we have 


n 2 


A + (77 + 2 ni) 2 I 


l>2n — Q {b>2n+2 + ^2n-2) — 0. 


(80) 
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Figure 10: (Color online) Diagram of amplitude versus frequency where hatched stable 
regions are found for a condensate containing a triply-charged vortex subjected to height 
modulation. We used Na s /d z = 125. Convergence is obtained already with two iterations. 
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At this point it is wise to define ladder operators L^ n 6 2 n = & 2 n ±2 which yields 

4 = |^ +jll + 2i (n ± l)] 2 1 - QIf„ ±2 | Q. 


(81) 


By using (81) to write (80) in terms of b 0 only, we obtain an iteration algorithm wherein we 


replace the ladder operator over and over inside itself which then becomes 

i -i 


A + ^if I - Q 


A+~^(v + 2*) 2 I 


+ 


A + 5-(r? - 2f) 2 / 


-i 


Q \b 0 = 0 . 


(82) 


Since we are not interested in trivial solutions for bo, the determinant of (82) must vanish. 


Thus the stability diagram for a modulation of the Gaussian potential with frequency D and 
amplitude Vo is presented in fig. [lOj where its resonant behavior does not depend on the initial 
conditions [3]. 
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The edges between stable and unstable domains (also called as Floquet fringes) were 
calculated by making rj — 0. Since the equilibrium configuration rarely has solution for 
Vo/hujp > Na s /d z , we only build the stability diagram for Vo/f\w p < Na s /d z . The iterative 
algorithm converges very fast, and does not require more than two iterations. 

The stable regions, also called resonance region, can lead the system to lose coherence if the 
excitation time is long enough (hundreds of milliseconds according to number of atoms) which 
leads to destruction of the condensate state. 

The dynamical mechanism works exciting the resonant mode by the oscillatory potential 
placed at the center of the condensate that suppresses completely the Qu-mode, when the 
correct frequency and amplitude are considered. Since this mode no longer exists, the vortex 


becomes stable (Fig. 11). It is what happens for the case where the static potential cannot 
stabilize the vortex by itself. On the other hand, in the case of static potential is enough to 
prevent the vortex decay, the modulation of the height plays an opposite role inducing the 
vortex decay in resonance regions. 


VIII. CONCLUSIONS 

In this paper we have studied the stability of collective modes as well as its dynamical 
stability for a quasi-2D Bose-Einstein condensate with a multi-charged vortex. The presence 
of a ^-charged vortex causes a shift in the frequencies of the cloud collective modes, however 
such changes are not substantial. The vortex rotational mode is an independent degree of 
freedom and does not affect vortex stability. The vortex dynamics couples with collective 
excitations, and it can be the cause for the ^-charged vortex decay. Its decay has as responsible 
the quadrupole oscillation Q v , which is one channel that leads the t-charged vortex to decay 
into £ singly vortices. This quadrupole is the main channel to doubly-charged vortex decay 
into two singly vortices. By applying a static Gaussian potential we can prevent the decay of 
a vortex for specific potential amplitudes, whereas for some regions in the parameter space can 
be stabilized by a time periodic modulation of the laser potential. 
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Figure 11: (Color online) Time evolution of condensate density with a triply-charged vortex 
for both free Gaussian potential (a, c, e) and dynamical potential (b, d, f). We used 
IJ,/huj p = 20, Na s /d z = 125, V 0 /hu p = 50, 5V/V 0 = 0.5, Vt/uj p = 5.2, e = 0.001, and a factor of 

0.01 multiplying the amplitude of deviations. 

Appendix A: Functions Aj (f, ao) and /*(£, «o) 

Similar functions to Aj (£, «o) for a 3D case have been calculated in Ref.(22J. Since it is a 
Thomas-Fermi wave-function the procedure to evaluate each integral is the same, where we 
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start changing the scale of both x and y coordinates according to x —>■ R x x and y —y R y y. By 
doing this £* becomes cp = £i/Ri, i.e. the integral becomes dimensionless. Now it is convenient 
to change the coordinates from cartesians to polar (x = p cos 0 and y = p sin &) where the 
integration domains are 0 < p < 1 and 0 < (j) < 2 tt, in this way we have 
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Where p F g (ai,..., a p ; b \,..., b q ] x ) are the hypergeometric functions, and B (x; a, 6) are beta 
functions. The functions derived from Gaussian potential have not an easy general form, then 
we write them in integral form: 
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